
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: SI6_additional_results.R
# Purpose: Produces results presented in section SI5
# Date: June 2025
# Data: pulled through 00_data_prep.R
#
# See 00_data_prep.R for technical disclaimer on R versions
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Load data
source("00_data_prep.R")


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Subset: Urban and rural respondents ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

table(full$urban, useNA="ifany")

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Urban respondents (Figure A9) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$urban==1,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
       hypothesis=="(1 Wales) - (1 Northern England)" |   
       hypothesis=="(2 Scotland) - (2 Northern England)" |
       hypothesis=="(2 Wales) - (2 Northern England)" |
       hypothesis=="(3 Scotland) - (3 Northern England)" |
       hypothesis=="(3 Wales) - (3 Northern England)" |
       hypothesis=="(4 Scotland) - (4 Northern England)" |
       hypothesis=="(4 Wales) - (4 Northern England)") %>%
mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
    grepl("Wales",hypothesis) ~ 4),
    group=case_when(grepl("1",hypothesis) ~ 1,
    grepl("2",hypothesis) ~ 2,
    grepl("3",hypothesis) ~ 3,
    grepl("4",hypothesis) ~ 4),
    p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
    p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
select(rowid,group,hypothesis,p5,p10) %>%
right_join(d2,by = join_by("rowid","group")) %>%
arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                         color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_urban.pdf", width=6, height=4)


# Rural respondents (Figure A10) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$urban==0,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_rural.pdf", width=6, height=4)



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Subset: Proximity to fossil fuel industry (Figure A11) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$ff_vicinity==1,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_ff_proximity.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (C) Subset: Climate change concern (Figure A12) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$cc_concern>1,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_cc_concern.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (D) Subset: Mitigation and adaptation ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

table(full$mitig_v_adapt, useNA="ifany")

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Mitigation (Figure A13) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$mitig_v_adapt==1 & is.na(full$mitig_v_adapt)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_mitigation.pdf", width=6, height=4)


# Adaptation (Figure A14) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$mitig_v_adapt==2 & is.na(full$mitig_v_adapt)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_adaptation.pdf", width=6, height=4)


# Both (Figure A15) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$mitig_v_adapt==3 & is.na(full$mitig_v_adapt)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_both.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (E) Subset: Economic security ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Below median income (Figure A16) ----
# Estimate model
table(full$income, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$income<=3 & is.na(full$income)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_below_median.pdf", width=6, height=4)



# Above median income (Figure A17) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$income>3 & is.na(full$income)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_above_median.pdf", width=6, height=4)


# Worry about local economy (Figure A18) ----
# Estimate model
table(full$worry_econdev, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$worry_econdev>=3 & is.na(full$worry_econdev)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_econ_worry.pdf", width=6, height=4)



# Vulnerable sectors (Figure A19) ----
# Estimate model
table(full$empl_sector, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$empl_sector<=6 & is.na(full$empl_sector)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_vulnerable_sector.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (F) Subset: Brexit vote ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Brexit vote as control (Figure A20) ----
# Create new leave variable
table(full$Brexit, useNA="ifany")
full <- full %>% 
  mutate(leave=case_when(Brexit=="leave" ~ 1, Brexit=="remain" ~ 0, .default = NA))

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)", "factor(leave)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full)

# Compare predicted probabilities for remainers and leavers
predictions(m2, newdata=datagrid(leave=c(0,1), sample=c("General population", "Northern England", "Scotland", "Wales"))) %>%
  group_by(group,sample) %>%
  select(rowid, group, sample, estimate, leave) %>%
  arrange(group, sample) %>%
  summarize(diff=estimate[leave==1]-estimate[leave==0])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_brexit_control.pdf", width=6, height=4)


# Remainers (Figure A21) ----
# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$leave==0,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_remainers.pdf", width=6, height=4)


# Leavers (Figure A22) ----

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$leave==1,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_leavers.pdf", width=6, height=4)

# Leavers and non-voters in 2019 GE (Figure A23) ----

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$leave==1 | full$voted=="no",])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_leavers_nonvoters.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (G) Subset: Devolution preferences ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Devolution preference (Figure A24) ----
# Estimate model
table(full$citizen, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[(full$citizen=="devolved" | full$citizen=="local") & is.na(full$citizen)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_devolution_pref.pdf", width=6, height=4)


# UK preference (Figure A25) ----
# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$citizen=="UK" & is.na(full$citizen)==FALSE,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_uk_pref.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (H) Subset: Vote choice ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Tory vote in 2019 election (Figure A26) ----
# Estimate model
table(full$parties, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$parties==1,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_Conservative_vote.pdf", width=6, height=4)


# Labour vote in 2019 election (Figure A27) ----
# Estimate model
table(full$parties, useNA="ifany")
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$parties==2,])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p10)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p10)==FALSE & p10=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI6_Labour_vote.pdf", width=6, height=4)



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                         END OF FILE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


